#include "SinglePhasePyrolysisMaterial.h"
#include<iostream>
#include "libmesh/quadrature.h"
using std::cout;
using std::endl;
registerMooseObject("TrilobitaApp", SinglePhasePyrolysisMaterial);

template<>
InputParameters validParams<SinglePhasePyrolysisMaterial>()
{
  InputParameters params = validParams<Material>();
  params.addParam<std::string>("rho_name", "Material property density name");
  //add regular thermal properties
  params.addRequiredParam<Real>("rhov_i","density of original material");
  params.addRequiredParam<Real>("rhoc_i","density of fully charred material");
  //add special parameters to compute pyrolysis Arrhenius relation
  params.addRequiredParam<Real>("T_begin", "temperature at which the reaction begins");
  params.addRequiredParam<Real>("precofficient","pre-exponential factor");
  params.addRequiredParam<Real>("ER","E/R");
  params.addRequiredParam<Real>("reaction_level","Reaction level of pyrolysis,normally 1,2 or 3");
  //add coupled variables
  params.addRequiredCoupledVar("temperature", "Coupled Temperature");
  params.addRequiredCoupledVar("rhoi", "Coupled density");
  return params;
}

Real SinglePhasePyrolysisMaterial::Source(Real rho, Real T, Real T_begin, Real rhoc, Real rhov, Real precoff, Real ER, Real reaction_level)
	{
	if (T>T_begin)
	{
		if ((rho-rhoc)/(rhov-rhoc)<0.01)
		{
			return 0;
		}
		else
		{
			Real tmp = (rho-rhoc)/rhov;
			return precoff*exp(-ER/T)*rhov*pow(tmp,reaction_level);
		}
	}
		else
	{
         return 0 ;
	}
}
SinglePhasePyrolysisMaterial::SinglePhasePyrolysisMaterial(const InputParameters & parameters) :
    	    Material(parameters),
		_rho_name(isParamValid("rho_name") ? getParam<std::string>("rho_name") : ""),
		_T_begin(declareProperty<Real>("T_begin"+_rho_name)),
		_precofficient(declareProperty<Real>("precofficient"+_rho_name)),
		_ER(declareProperty<Real>("_ER"+_rho_name)),
		_reaction_level(declareProperty<Real>("reaction_level"+_rho_name)),
		_rhov_i(declareProperty<Real>("rhov_i"+_rho_name)),
		_rhoc_i(declareProperty<Real>("rhoc_i"+_rho_name)),
		_rho_i(declareProperty<Real>("rho_i"+_rho_name)),
		_pyrolysis_source(declareProperty<Real>("pyrolysis_source"+_rho_name)),
		_dsource_drho(declareProperty<Real>("dsource_drho"+_rho_name)),
		_dsource_dT(declareProperty<Real>("dsource_dT"+_rho_name)),
		_T_value(coupledValue("temperature")),
		_gradient_T_value(coupledGradient("temperature")),
		_rho_value(coupledValue("rhoi")),
		_rho_dt_value(coupledDot("rhoi")),
		_rhodot_drho_value(coupledDotDu("rhoi"))
{
	_rho_i_v_value =  getParam<Real>("rhov_i");
	_rho_i_c_value =  getParam<Real>("rhoc_i");
	_T_begin_value =  getParam<Real>("T_begin");
	_precoff_value =  getParam<Real>("precofficient");
	_ER_value =  getParam<Real>("ER");
	_reaction_level_value =  getParam<Real>("reaction_level");
}

void SinglePhasePyrolysisMaterial::computeProperties()
{
	for (unsigned int qp(0);qp < _qrule->n_points(); ++qp)
	{
//compute regular thermal properties for orinin and charred material separately
	_rhov_i[qp] = _rho_i_v_value;
	_rhoc_i[qp] = _rho_i_c_value;
//compute special parameters to caculate pyrolysis Arrhenius relation
	_T_begin[qp] = _T_begin_value;
	_precofficient[qp]  = _precoff_value;
	_ER[qp]  = _ER_value;
	_reaction_level[qp]  = _reaction_level_value;
//compute solid density material parameter
	_rho_i[qp]  = _rho_value[qp];
//compute source
	Real epsi = 1E-8;
	Real source  = Source(_rho_value[qp], _T_value[qp],  _T_begin_value,  _rho_i_c_value,  _rho_i_v_value,  _precoff_value,  _ER_value, _reaction_level_value);
	Real source_new1 = Source(_rho_value[qp], _T_value[qp]+epsi,  _T_begin_value,  _rho_i_c_value,  _rho_i_v_value,  _precoff_value,  _ER_value, _reaction_level_value);
	Real source_new2 = Source(_rho_value[qp]+epsi, _T_value[qp],  _T_begin_value,  _rho_i_c_value,  _rho_i_v_value,  _precoff_value,  _ER_value, _reaction_level_value);
	_pyrolysis_source[qp] = Source(_rho_value[qp], _T_value[qp],  _T_begin_value,  _rho_i_c_value,  _rho_i_v_value,  _precoff_value,  _ER_value, _reaction_level_value);

	_dsource_dT[qp] = (source_new1-source)/epsi;
	_dsource_drho[qp] = (source_new2-source)/epsi;

	}

}


